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Determination of ^ ^ 0.8 neutral hydrogen fluctuations using the 
21 cm intensity mapping auto-correlation 
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ABSTRACT 

The large-scale distribution of neutral hydrogen in the universe will be luminous through its 
21 cm emission. Here, for the first time, we use the auto-power spectrum of 21 cm intensity 
fluctuations to constrain neutral hydrogen fluctuations at z ^ 0.8. Our data were acquired 
with the Green Bank Telescope and span the redshift range 0.6 < z < 1 over two fields total- 
ing « 41 deg. sq. and 190 hr of radio integration time. The dominant synchrotron foregrounds 
exceed the signal by ~ lO'^, but have fewer degrees of freedom and can be removed effi- 
ciently. Even in the presence of residual foregrounds, the auto-power can still be interpreted 
as an upper bound on the 2 1 cm signal. Our previous measurements of the cross-correlation of 
21 cm intensity and the WiggleZ galaxy survey provide a lower bound. Through a Bayesian 
treatment of signal and foregrounds, we can combine both fields in auto- and cross-power 
into a measurement of r^Hi^Hi = [0.62^Q^g] x 10^'^ at 68% confidence with 9% systematic 
calibration uncertainty, where Ohi is the neutral hydrogen (H l) fraction and &hi is the H I 
bias parameter We describe observational challenges with the present dataset and plans to 
overcome them. 
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1 INTRODUCTION 

There is substantial interest in the viability of cosmological struc- 
ture surveys that map the intensity of 21 cm emission from neu- 
tral hydrogen. Such surveys could be used to study large-scale 
structure (LSS) at intermediate redshifts, or to study the epoch 
of reionization at high redshift. Surveys of 21 cm intensity have 
the potential to be very efficient since the resolution of the in- 
strume nt can be matched to the large scales of cosmological in- 
terest dCh ang et al.'l 2008l ; iLoeb & Wvith3 l2008l : ISeo et al.ll201(]| : 
lAnsari et a l. 2012a). Several exp eriments, including BAOBAB 



Poberetal.ll 2013b). BAORadio JAnsari et al.l l2012bl) BINGO 
Battve et aljl2012i) . CHIM^E and TianLai ( ICherj|2012l) propose 
to conduct redshift surveys from z ~ 0.5 to z ~ 2.5 using this 



method. 



* E-mail: eswitzer@cita.utoronto.ca 
f E-mail: kiyoScita.utoronto . ca 
|http: //chime .phas .ubc ■ ca"71 



The principal challenges for 21 cm experiments are astronom- 
ical foregrounds and terrestrial radio frequency interference (RFI). 
Extragalactic sources and the Milky Way produce synchrotron 
emission that is three orders of magnitude brighter than the 21 cm 
signal. However, the physical process of synchrotron emission is 
known to produce spectrally-smooth radiation, occupying few de- 
grees of freedom along each line of sight. In the absence of instru- 
mental effects, these degrees of freedom are thought to be separa- 
ble from the signal (JLiu & Tegmark 201 1, 2012; Shaw et al. 2013]). 
RFI can be minimized through site location, sidelobe control, and 
band selection. In the Green Bank Telescope data analyzed here, 
RFI is not found to be a significant challenge or limiting factor. 

Subtraction of synchrotron emission has proven to be chal- 
lenging in practice. Instrumental effects such as passband calibra- 
tion and polarization leakage couple bright foregrounds into new 
degrees of freedom that need to be removed from each line of sight 
to reach the level of the 21cm signal. The spectral functions de- 
scribing these systematics can not all be modeled in advance, so 
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we take an empirical approach! to foreground removal by estimat- 
ing dominant modes from the covariance of the map itself. This 
method requires more caution because it also removes cosmologi- 
cal signal, which must be accounted for. 

Large-scale neutral hydrogen fluctuations above redshift 
z = 0.1 have been unambiguously detected only in cross- 
correlatiori with existing siirveys of optically-selec ted galaxies 
JLah et al.l2009l : IChang et alj2010l . lMasui et alhoiSh . Here, resid- 
ual 21 cm foregrounds boost the errors but do not correlate with the 
optical galaxies. The density fluctuations traced by survey galaxies 
may not coirelate perfectly with the emission of neutral hydrogen, 
so their cross-correlation can be interpreted as a lower limit on the 
fluctuation power of 21 cm emission. 

Several efforts have used the 21 cm line to pla c e up- 
per bounds on the reionization era (JBebbingtorJ 1986t 
iBowman & Rogers! |2010| : Paciga e t'ai]|2013l; [P ober et aj 1201 3al) 
and z ^ 3 (see e.g., IWieringa. de Bruvn & Katgerli ( Il992f) : 
ISubrahmanvan & Anantharamaiahl 1 119901) ^ without the need to 
cross-correlate with an external data set. This is the first work to 
describe similar bounds for z ~ 0.8, using two fields totaling 
« 41 deg. sq. and 190 hr of radio integration time with the Green 
Bank Telescope. Unlike the bounds from reionization, for which 
there is currently no cross-correlation, we are able to combine 
the auto- and cross-powers in a novel way, making a Bayesian 
inference of the amplitude of neutral hydrogen fluctuations, 
parameterized by fiHibHi- 

Throughout, w e use cosmological parameters from 
iKomatsu etal] J2009I') . 



2 OBSERVATIONS AND ANALYSIS 

The analysis here is based on th e same observation s used for the 
cross-correlation measurement in lMasui et al.l ( 120131) . We flag RFI 
in the data, calculate 3D intensity map volumes, clean foreground 
contamination, and estimate the power spectrum. Here we will 
summarize essenti al aspects of the observations and analysis in 
iMasui et al.l ( l2013r) . and describe the auto-power analysis in more 
detail. 

Observations were conducted with the 680 — 920 MHz prime- 
focus receiver at the Green Bank Telescope (GBT), sampled from 
700 MHz (2 = 1) to 900 MHz {z = 0.58) in 256 uniform spectral 
bins. The analysis here uses a 105 hr integration of a 4.5° x 2.4° 
15 hr "deep" field centered on 14''31™28.5= right ascension, 2°0' 
declination and an 84 hr integration on a 7.0° x 4.3° Ihr "wide" 
field centered on 0'^52™0^ right ascension, 2°9' declination. 

The beam FWHM at 700 MHz is 0.314° , and at 900 MHz it is 
0.250° . At band-center, the beam width corresponds to a comoving 
length of 9.6ft^^Mpc. Both fields have nearly complete angular 
overlap and good redshift cover age with the WiggleZ Dark Energy 
Survey JDrinkwater et al.ll2010l) . Our absolute calibrat ion is deter- 
mined from radio point sources and is accurate to 9% l IMasui et alj 
l2013h . For clarity, this remains as a separately quoted systematic 
error throughout, and plotted posterior distributions are based on 
statistical eiTors only. 



2.1 Foreground Cleaning 

In this section, we develop the map cleaning formalism and dis- 
cuss its connection to survey strategy. Begin by packing the three- 
dimensional map into an A^^ x Ng matrix M by unwrapping the 
Ng RA, Dec pointings. For the moment, ignore thermal noise 



in the map. The empirical v — v' covariance of the map is 
C = MM^^/TVe, and it contains both foregrounds and 21 cm sig- 
nal. This can be factored as C = UAU"'^, where A is diagonal and 
sorted in descending value. From each line of sight, we can then 
subtract a subset of the modes U that describe the largest compo- 
nents of the variance through the operation (1 — USU^)M, where 
S is a selection matrix with 1 along the diagonal for modes to be 
removed and elsewhere. 

In reality, M also contains thermal noise. To minimize its in- 
fluence on our foreground mode determination, we find the noise- 
inverse weighted cross-variance of two submaps from the full sea- 
son of observing. Here, Cab = (Wa o Ma)(Ws o Ms)^/A''e, 
where A and B denote sub-season maps, Wa is the noise inverse- 
variance weight per pixel of map A (neglecting correlations), and 
o is the element-wise matrix product. Cab is no longer symmet- 
ric, and we take its singular value decomposition (SVD) instead, 
using the left and right singular vectors to clean maps A and B re- 
spectively. The weights are calculated in the noise model developed 
in the map-maker, but roughly track the map's integration depth 
and weigh against RFI. The weight is nearly separable into angle 
(through integration time) and frequency (through Tsysiy)), but we 
average to make it formally separable and so rank-1, so that it does 
not increase the map rank. The weighted removal for map A be- 
comes (1/Wa) o (1 - UaSU3J)Wa o Ma, where 1/Wa is the 
element-wise reciprocal. 

Our empirical approach to foreground removal is limited by 
the amount of information in the maps. The fundamental limitation 
here surprisingly is not from the number of degrees of freedom 
along the line of sight, but is instead t he number of inde pendent 
angular resolution elements in the map (lNitvanandall20 1 Of) . To see 
why this is the case, notice that in the absence of noise, our clean- 
ing algorithm is equivalent to taking the SVD of the map directly: 
M = USV^ and thus C oc MM"^ = US^U"^, with the same 
set of frequency modes U appearing in both decompositions. The 
rank of C coincides with the rank of M and is limited by the num- 
ber of either angular or frequency degrees of freedom. 

Taking the foreground modes to all have comparable spurious 
overlap with the signal, one anives at a transfer function rule of 
thumb T = PBig.out/Psig.in ~ (1 - Ar,„/iV,)(l - Ar™/iV,es), 
where 7V,„ is the number of modes removed, N^ = 256 is the 
number of frequency channels and A'^res is the number of resolu- 
tion elements (roughly the survey area divided by the beam solid 
angle). A limited number of resolution elements can greatly reduce 
the efficacy of the foreground cleaning at the expense of signal. 

The wide and deep fields respectively have central low-noise 
regions of ~ 10 deg. sq. and ~ 3 deg. sq., giving roughly 90 and 
30 independent resolution elements at the largest beam size. The 
rank of C is then less than the number of available spectral bins in 
both cases. To increase the effective number of angular degrees of 
freedom that enter the weighted v — v' covariance, we saturate the 
angular weight maps at their 50'th percentile. 

The choice of the number of modes to remove is a trade-off be- 
tween bias from residual foregrounds (for too few modes removed), 
and increasing errors (from signal lost as too many modes are re- 
moved). We find that 30 (10) modes for the wide (deep) field shows 
a good balance of minimal foregrounds and errors. 

2.2 Instrumental Systematics 

The physical mechanism of synchrotron radiation suggests that it 
is des cribed by a handful of smooth modes along each line of 
sight JLiu & Tegmarkll2012h . Instrumental response to bright fore- 
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grounds, however, can convert these into new degrees of free- 
dom. An imperfect and time-dependent passband calibration will 
cause intrinsically spectrally smooth foregrounds to occupy multi- 
ple modes in our maps with non-trivial spectral structure. We con- 
trol this using a pulsed electronic calibrator, averaged for each scan. 

We believe that the most pernicious spectral structure is 
caused by leakage of polarization into intensity. Our observed po- 
larization mixing is ~ 10% between Stokes parameters and has 
both angular and spectral structure. The spectral structure converts 
spectrally smooth polarization into new degrees of freedom. Fara- 
day rotation of the polarization introduces further spectral degrees 
of freedom. The majority of the polarization leakage is observed 
to be an odd function on the sky, slightly broader than the primary 
intensity response beam. To mitigate this leakage, we convolve the 
maps to a common resolution corresponding to 1.4 times the beam 
size at 700 MHz (the largest beam). This also treats the leakage 
of spatial structure into spectral structure from frequency depen- 
dence of the beam. Such a convolution is viable because GBT has 
roughly twice the resolution needed to map large-scale structure in 
the linear regime. However, this convolution reduces the number 
of independent resolution elements in the map by a factor of two, 
increasing the challenges discussed in Sec. 12. II 

The present results are limited largely by the area of the re- 
gions and our understanding of the instrument. With a factor of 
roughly ten more area, the resolution could be degraded at less 
expense to the signal. This requires significant telescope time be- 
cause the area must also be covered to roughly the same depth as 
our present fields. It would however provide a significant boost in 
overall sensitivity for scientific goals such as measurement of the 
redshift-space distortions. In addition, we are investigating map- 
making that would unmix polarization using the Mueller matrix as 
a function of offset from the boresight, as determined from source 
scans. 



2.3 Power Spectrum Estimation 

Our starting point for power spectr al estimation is the opti- 
mal quadratic estimator described in l ILiu & TegmarkI EOl lb . To 
avoid the thermal noise bias , we only consider c ross-powers be- 
tween four sub-season maps JTristram et al.ll2005h . labeled here as 
{A,B,C,D}. Thermal noise is uncorrelated between these sec- 
tions, which we have chosen to have similar integration depth and 
coverage. The foreground modes are determined separately for 
each side of the pair using the SVD of Sec. 12.11 Up to a normal- 
ization, the resulting estimator for the pair of submaps A and B is 
P{}^i)AxB oc (wAnAmA)^QiWsnBms. Here, we have un- 
wrapped the map matrix Ma into a one-dimensional map vector 
niA and written the foreground cleaning projection (1/Wa) o (1 — 
UaSUa) Wa o Ma as YIa'co^a- The weighted mean of each fre- 
quency slice of the map is also subtracted. The map weight wa 
is the matrix Wa used in the SVD, but unwrapped, and along the 
diagonal. Procedurally, the estimator amounts to weighting both 
foreground-cleaned maps, taking the Fourier transform, and then 
summing the three-dimensional cross-pairs to find power in annuli 
in two-dimensional k space, k; = {fcx,i, fc[|,i}- The Fourier trans- 
form and binning are performed by Qi here. We calculate six such 
crossed pairs from the four-way subseason split of the data, and let 
the average over these be the estimated power -P(ki). 

We calculate transfer functions to describe signal lost in the 
foreground cleaning and through the finite instrumental resolution. 
These are functions of fcx and fcy. The beam transfer function is 
estimated using Gaussian 21 cm signal simulations that have been 



convolved by the beam model. The foreground cleaning transfer 
function can be efficiently estimated through Monte Carlo simula- 
tions as 



T(kO = 



[vyAnA+s(mA + ms) - w AnAVHA]^ Qi^m^ 

(wrAms)^Q,ms 



,(1) 



where the A-\-s subscript denotes the fact that the foreground clean- 
ing modes have been estimated from SiV — v' covariance that has 
added 21 cm simulation signal, nis. The limited number of angu- 
lar resolution elements (Sec. 12.1) results in an anticorrelation of 
the cleaned foregrounds with the signal itself, represented by the 
term (wAnA+smA)"^Qinis. To reduce the noise of the simula- 
tion cross-power, note that we subtract ■waHavha in the numera- 
tor. Finally, we find the weighted average of these across the four- 
way split of maps. We find that 300 signal simulations are sufficient 
to estimate the transfer function. 

After compensating for lost signal using transfer functions for 
the beam and foreground cleaning, we bin the two-dimensional 
powers onto one-dimensional band-powers. We weight bins 
by their two-dimensional Gaussian inverse noise variance ex 
iV(ki)T(ki)^/Pauto(ki)^, where Pauto(ki) is the average of 
{PaxA,Pbxb,PcxC,Pdxd} (pairs which contain the thermal 
noise bias), and N(\i.i) is the number of three-dimensional k cells 
that enter a two-dimensional bin ki. In addition to the Gaus- 
sian noise weights, we impose two additional cuts in the two- 
dimensional k power. For fc|| < 0.035 /i/Mpc, k± < 0.08 /i/Mpc 
for the deep field, and fcx < 0.04 fe/Mpc for the wide field, there 
are few harmonics in the volume, resulting in "straggler" strips in 
the two-dimensional power spectrum where the errors are poorly 
estimated. For fcx > 0.3 h/Mpc, the instrumental resolution pro- 
duces significant signal loss, so this is truncated also. 

Foregrounds in the input maps and the 21 cm signal itself are 
non-Gaussian, but after cleaning, the thermal noise dominates both 
contributions in an individual map, and Gaussian errors (see, e.g. 
lOas et alj ( 1201 ih ) provide a reasonable approximation. These take 
as input the auto-power measurement itself (for sample variance) 
and PaxA terms that represent the thermal noise. Sample variance 
is significant only in the deep field in the lower 1/3 of the reported 
wavenumbers. Gaussian eiTors agree with the standard deviation of 
the six crossed-pairs that enter the spectral estimation in the regime 
where sample variance is negligible. 

The finite survey size and weights result in correlations be- 
tween adjacent k bins. We apodize in the frequency direction using 
a Blackman window, and in the angular direction using the map 
weight itself (which falls off at the edges due to scan coverage). 
The bin-bin correlations are estimated using 3000 signal plus ther- 
mal noise simulations assuming T^ys = 25 K. To construct a full 
covariance model, these are then re-calibrated by the outer product 
of the Gaussian error amplitudes for the data relative to the thermal 
noise simulation errors. 

The Bayesian method developed in the next section assumes 
that adjacent bins are uncorrelated. To achieve this, we take the 
matrix square-root of the inverse of our covariance model ma- 
trix and normalize its rows t o sum to one. This provides a set of 
functions which decorrelates ( [Hamilton & TegmarkiEOOOr) the pre- 
whitened power spectrum and boosts the errors. At large scales 
(k — 0.1 h/Mpc) where these effects are relevant, decorrelation 
and sample variance increase the errors by a factor of 1.5 in the 
wide field and 4 in the deep field. 

The methods that we have developed for calculating maps 
from data and power spectra from maps will be described and more 
fully motivated in a future paper. 
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Figure 1. Temperature scales in our 21 cm intensity mapping survey. The 
top curve is the power spectram of the input deep field with no cleaning ap- 
plied (the wide field is similar). Throughout, the deep field results are green 
and the wide field results are blue. The dotted and dash-dotted lines show 
theiTnal noise in the maps. The power spectra avoid noise bias by crossing 
two maps made with separate datasets. Nevertheless, thermal noise limits 
the fidelity with which the foreground modes can be estimated and removed. 
The points below show the power spectram of the deep and wide fields af- 
ter the foreground cleaning described in Sec. 12. II Negative values are shown 
with thin lines and hollow markers. Any residual foregrounds will additively 
bias the auto-power The red dashed line shows the 21cm signal expected 
from the amplitude of the cross-power with the WiggleZ survey (for r = 1) 
and based on simulations processed by the same pipeline. 



3 RESULTS 

The auto-power spectra presented in Figure [T] will be biased by 
an unknown positive amplitude from residual foreground contam- 
ination. These data can then be interpreted as an upper bound 
on the neutral hydrogen fluctuation amplitude, fiHi^Hi- In addi- 
tion, we have also m easured the cross-c orrelation with the Wig- 
gleZ Galaxy Survey JMasui et al.ll2013f) . This finds Q,uibuir = 
[0.43 ± 0.07(stat.) ± 0.04(sys.)] x 10"^, where r is the Wig- 
gleZ galaxy-neutral hydrogen cross-correlation coefficient (taken 
here to be independent of scale). Since \r\ < 1 by definition and is 
measured to be positive, the cross-correlation can be interpreted as 
a lower bound on Quibni- In this section, we will develop a pos- 
terior distribution for the 2 1 cm signal auto-power between these 
two bounds, as a function of k. We will then combine these into a 
posterior distribution on Quibm. 

The probability of our measurements given the 21 cm signal 
auto-power and foreground model parameters is 

pid,\e,)=p{d^\s,,r)pi4^^^\s,,ft^nPidT'''\sk,fr''l. (2) 

Here, dfc — {d'^, dj,'^°'', d""^°} contains our cross-power and 
deep and wide field auto-power measurements, while Ok = 
{sk, r, fk"'''^, /ir"*"^} contains the 21 cm signal auto-power, cross- 
coiTelation coefficient, and deep and wide field foreground con- 
tamination powers, respectively. The cross-power variable d'^ rep- 
resents the constraint on Quibnir from b oth fields and the range of 
wavenumbers used in lMasui et alj ( 120130 . The band-powers d^.*^^^ 
and dj^"^^ are independently distributed following decorrelation of 
finite-survey effects. We assume that the foregrounds are uncorre- 
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Figure 2. Comparison with the thermal noise limit. The dark and light 
shaded regions are the 68% and 95% confidence intervals of the measured 
21 cm fluctuation power The dashed fine shows the expected 21 cm signal 
implied by the WiggleZ cross-correlation if r = 1. The solid line represents 
the best upper 95% confidence level we could achieve given our eiTor bars, 
in the absence of foreground contamination. Note that the auto-correlation 
measurements, which constrain the signal from above, are unc orrelated be- 
tween k bins, while a single global fit to the cross-power (in iMasui et al J 
(20131)) is used to constrain the signal from below. Confidence intervals 
do not include the systematic caUbration uncertainty, which is 18% in this 
space. 



lated between k bins and fields, also. This is conservative because 
knowledge of foreground coiTelations would yield a tighter con- 
straint. We take p{d'^\sk, r) to be normally distributed with mean 
proportional to Vy/Ik, and p(d^°'^''|sfc, fjJ^'^'^) to be normally dis- 
tributed with mean Sk + f^^'^'^ ^nd errors determined in Sec 12.31 
(and analogously for the wide field). Only the statistical uncertainty 
is included in the width of the distributions, as the systematic cali- 
bration uncertainty is perfectly correlated between cross- and auto- 
power measurements and can be applied at the end of the analysis. 
We apply Bayes' Theorem to obtain the pos- 
terior distribution for the parameters, p{Ok\Ak) oc 
p(dfc|0fc)p(sfe)p(r-)p(/^™'')p(/""^°). For the nuisance pa- 
rameters, we adopt conservative priors. p{fiJ"^^) ^nd p(/fc'"^°) 
are taken to be flat over the range < /fc < oo. Likewise, we 
take p{r) to be constant over the range < r < 1, which is 
conservative given the theoretical bias toward r ~ 1. Our goal is 
to marginalize over these nuisance parameters to determine Sfe. We 
choose the prior on Sk, p{sk), to be flat, which translates into a 
prior p(S1hi&hi) oc Quibui- The data likelihood adds significant 
information, so the outcome is robust to choices for the signal 
prior. The signal posterior is 



p(sfc|dfc 



\ / / /"deep /■wide I 1 

) = / P(Sk,r,fk "^Jk \dk 



, 1 n /■deep 1 rwide 



(3) 



This involves integrals of the form J^ p(d''\s,r)p{r)dr which, 
given the flat priors that we have adopted, can generally be writ- 
ten in terms of the cumulative distribution function of p{d'^\s, r). 
Figure|2]shows the allowed signal in each spectral fc-bin. 

Taking the analysis further, we com bine band-powers i nto a 
single constraint on QnibHi- Following iMasui et alj ( l2013r) . we 
consider a conservative k range where errors are better estimated 
(fc > 0.12 /i/Mpc, to avoid edge effects in the decorrelation op- 
eration) and before uncertainties in nonlinear structure formation 
become significant (k < 0.3 /i/Mpc). Figure[3]shows the resulting 
posterior distribution. 

Our analysis yields ^mbm = [0.62l°;^^] x 10"^ at 68% 
confidence with 9% systematic calibration uncertainty. Note that 
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Figure 3. The posterior distribution for the parameter SIhi^hi coming from 
the WiggleZ cross-power spectrum, deep field and wide field auto-powers, 
as well as the joint likelihood from all three datasets. The individual distri- 
butions from the cross-power and auto-powers are dependent on the prior 
on nHl^Hl while the combined distribution is essentially insensitive. The 
distributions do not include the systematic calibration uncertainty of 9%. 



we are unable to calculate a goodness-of-fit to our model because 
each measurement is associated with a free foreground parameter 
which can absorb any anomalies. 



4 DISCUSSION AND CONCLUSIONS 

Through the measurement of the auto-power, we extend our pre- 
vious cross-power measurement of SIhi&hi'' JMasui et alj|2013l) 
to a determination of fiHi&Hi- This is the first constraint on the 
amplitude of 21 cm fluctuations at z ~ 0.8, and circumvents the 
degeneracy with the cross-correlation r. The 21cm auto-power 
yields a true upper bound because it derives from t he integral of 



the mass function. In the future, redshift distortions dWvithi 



;gral ot 
320081 : 



iMasui. McDonald & Penl2010l) can be used to further break the de- 
generacy between 6hi and Ohi, and complement cha llenging HST 
measurements of JIhi ( IRao. Tumshek & Nestoj200q) . Our present 
survey is limited by area and sensitivity, but we have shown that 
foregrounds can be suppressed sufficiently, to nearly the level of 
the 21 cm signal, using an empirical mode subtraction method. Fu- 
ture surveys exploiting the auto-power of 21 cm fluctuations must 
develop statistics that are robust to the additive bias of residual 
foregrounds, and control instrumental systematics such as polar- 
ized beam response and passband stability. 
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